### Latent Hawkishness Variable
install.packages("MCMCpack")
install.packages("doBy")
library(MCMCpack)
library(doBy)
set.seed(5)


## US Presidential Hawkishness
d <- read.csv("CarterSmithThetas-sample.csv", header=TRUE)
row.names(d) <- d$ctry_yr


d$theta2_mean <- ordered(as.factor(d$theta2_mean))
d$yr1953 <- ordered(as.factor(d$yr1953))
d$yr1961 <- ordered(as.factor(d$yr1961))
d$yr1963 <- ordered(as.factor(d$yr1963))
d$yr1969 <- ordered(as.factor(d$yr1969))
d$yr1974 <- ordered(as.factor(d$yr1974))
d$yr1977 <- ordered(as.factor(d$yr1977))
d$yr1981 <- ordered(as.factor(d$yr1981))
d$yr1989 <- ordered(as.factor(d$yr1989))
d$yr1993 <- ordered(as.factor(d$yr1993))
d$yr2001 <- ordered(as.factor(d$yr2001))
d$yr2005 <- ordered(as.factor(d$yr2005))

# Estimate model
post_us <- MCMCmixfactanal(~
                             theta2_mean+yr1953+yr1961+yr1963+yr1969
                           +yr1974+yr1977+yr1981+yr1989+yr1993++yr2001+yr2005,
                           factors=1, data=d,
                           lambda.constraints = list(theta2_mean=list(2,"+")),
                           burnin=20000, mcmc=100000, thin=50,
                           verbose=500, L0=.25, store.lambda=TRUE,
                           store.scores=TRUE, tune=1.2)

summary(post_us)

# convert object from mcmc to data.frame and transpose
c<-as.data.frame(post_us)
c<-t(c)
c1<-apply(c,1,quantile,probs = c(0.025, 0.05, 0.5, 0.95, 0.975))
mean<-apply(c,1,mean)
sd<-apply(c,1,sd)
c1<-rbind(c1,mean,sd)
c1<-t(c1)

# isolate and extract lambda/psi (beginning and end of variables)
lambda<-c1[1:44,]
n<-nrow(c1)
psi<-c1[(n-10):n,]
lambda<-rbind(lambda,psi)

# isolate country ideal points (phi; in middle)
phi<-c1[(n-10):n,]

# change vars to lower 95+ 90+ median+ and upper 90+ 95; change names to remove ".2"
lambda
colnames(lambda) <- c("low95","low90","median","high90","high95","mean","sd")

# change vars to lower 95+ median+ and upper 90; change names to leave just country name
phi
colnames(phi) <- c("low95","low90","median","high90","high95","mean","sd")
rownames(phi) <- sub("\\phi.", "", rownames(phi))

# Save lambda/psi as one .csv file and save country ideal points as another .csv file
write.csv(lambda,file="lambda_us.csv")
write.csv(phi,file="country_irt_us.csv")

